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Abstract 

When a system has more than one stable state, how can the stability of these states be 
compared? This deceptively simple question has important consequences for ecosystems, 
because systems with alternative stable states can undergo dramatic regime shifts. The 
probability, frequency, duration, and dynamics of these shifts will all depend on the relative 
stability of the stable states. Unfortunately, the concept of “stability” in ecology has 
suffered from substantial confusion and this is particularly problematic for systems where 
stochastic perturbations can cause shifts between coexisting alternative stable states. A 
useful way to visualize stable states in stochastic systems is with a ball-in-cup diagram, in 
which the state of the system is represented as the position of a ball rolling on a surface, 
and the random perturbations can push the ball from one basin of attraction to another. 
The surface is determined by a potential function, which provides a natural stability metric. 
However, systems amenable to this representation, called gradient systems, are quite rare. 
As a result, the potential function is not widely used and other approaches based on linear 
stability analysis have become standard. Linear stability analysis is designed for local 
analysis of deterministic systems and, as we show, can produce a highly misleading picture 
of how the system will behave under continual, stochastic perturbations. In this paper, we 
show how the potential function can be generalized so that it can be applied broadly, 
employing a concept from stochastic analysis called the quasi-potential. Using three classic 
ecological models, we demonstrate that the quasi-potential provides a useful way to 
quantify stability in stochastic systems. We show that the quasi-potential framework helps 
clarify long-standing confusion about stability in stochastic ecological systems, and we 
argue that ecologists should adopt it as a practical tool for analyzing these systems. 


Keywords: alternative stable states, stochastic dynamics, regime shifts, quasi-potential, 
Freidlin-Wentzcll, stochastic differential equations, Hamilton-Jacobi, resilience 
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Introduction 

Researchers have long been fascinated by the possibility for ecosystems to have more 


than one stable state (May 1977, Beisner et ah 2003). Such ecosystems have been observed 


in both natural (van de Koppel et ah 2001) and experimental (Chase 2003) settings. 
Systems with multiple (i.e., alternative) stable states can can abruptly shift from one 


stable state to another, sometimes with catastrophic consequences (Scheffer and Carpenter 


2003), so understanding their properties is crucially important. 


Unfortunately, the understanding of alternative stable states has been significantly 
hampered by ambiguity about the term “stable”. Grimm and Wissel (2008) note that 
stability is “one of the most nebulous terms in the whole of ecology,” and they catalog 163 
different definitions. Much of this confusion arises when researchers attempt to apply tools 
designed for the analysis of deterministic models to stochastic models. Fortunately, there is 


a well-developed mathematical framework, the Freidlin-Wentzell quasi-potential (Freidlin 


and Wentzcll 2012), that provides a rigorous yet natural way to understand alternative 


stable states in stochastic systems. In this paper, we explain how this tool can clarify much 
of the confusion about stability in ecological systems by translating intuitive concepts into 
quantifiable mathematical properties. Through three examples, we show how the quasi¬ 
potential serves as a useful metric of stability, and allows for effective stability comparison 
between alternative stable states. The results from quasi-potential analysis often contrast 
with those from standard stability analysis, and our examples explore these discrepancies. 
Furthermore, the quasi-potential allows for stability to be quantified on a continuum that 
corresponds well with the system’s dynamics, and it can be applied to any system state, 
regardless of whether that state is a deterministic equilibrium. Using the quasi-potential, a 
system can be decomposed into orthogonal components, and we explain how this 
decomposition can be interpreted ecologically. Finally, the quasi-potential offers insight 
into the most probable paths a system will take in transitioning from one state to another. 

Holling’s foundational work on resilience and stability anticipated the quasi-potential’s 
basic essence ( |Holling 1973); later, Tuljapurkar and Semura (1979) made the insight that 
Holling’s intuitive ideas were connected to the mathematical work of Freidlin and Wentzell 
(1970). At that time, numerical methods were insufficient to allow for general, practical 


computation of quasi-potentials (see Ludwig 1975), so Tuljapurkar and Semura’s insight 


did not receive the recognition it deserved. In subsequent decades, the flurry of research on 
alternative stable states largely overlooked this insight. Recently, the quasi-potential has 
been embraced by researchers analyzing models in other areas of biology, although it often 
appears under other names, and is disconnected from the Freidlin-Wentzell formulation 


(but see Zhou 2012). These applications include gene regulatory networks (Lv et al. 2014 


Zhou et al. 2012), neural networks (Yan et al. 2013), and evolution (Zhang et al. 2012 
Wang et al. 2011). Very recently, it has been applied to a predator-prey system (|Xu et al. 


2014), and with countless other possibilities for application, we argue that the quasi¬ 


potential is poised to become a major quantitative tool in ecology. 

This paper makes three novel contributions to the field of ecology. First, it shows how 
the quasi-potential can clarify the confusing tangle of stability concepts that confront 
ecologists. Second, it demonstrates how the quasi-potential can be used to quantify 
stability in systems with alternative stable states, and how the results can be different from 
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and often more useful than deterministic methods. Finally, it shows how a new numerical 


algorithm for the computation of quasi-potentials (Cameron 2012) can be expanded for 


application to systems with multiple stables states, and highlights the utility of the quasi¬ 
potential for understanding such systems. 

We use three well-established ecological models to illustrate these ideas. First, we show 
how traditional linear stability analysis fails to capture the salient features of a stochastic 
lake eutrophication model, and explain how the system’s potential function provides more 
useful analytic insights. Next, we move to higher-dimensional systems, where potential 
functions rarely exist. We explore a consumer-resource model with alternative stable states 
that does not have a potential function. We explain how the quasi-potential is defined, and 
show its usefulness in analyzing this model. Finally, we explore another consumer-resource 
model with a stable limit cycle to demonstrate how the quasi-potential is useful when stable 
states are more complicated than point equilibria. We conclude by discussing the quasi¬ 
potential as a unifying framework for existing notions of stability in stochastic systems. 

Example 1: Lake Eutrophication 

Lake ecosystems are among the most well-studied examples of alternative stable states 
in ecology. A foundational model by Carpenter et al. (1999) successfully describes the 
coexistence of a eutrophic state, corresponding to high phosphorous concentration, and an 
oligotrophic state, corresponding to low phosphorous concentration. Later work by Guttal 


and Jayaprakash (2007) showed how stochasticity can cause this system to switch between 


the two stable states, and we will use their model as a starting point for exploring the 
quantification of stochastic stability. 

The underlying deterministic model (i.e., the “deterministic skeleton”) describes how 
the nutrient (phosphorous) concentration x changes over time: 


dx 

dt 


c — sx + r 


x H 


Xq + x q 


( 1 ) 


c is the nutrient inflow rate and s is the nutrient loss rate (due to sedimentation, outflow, 
and sequestration in benthic plants). The last term represents nutrient recycling, r is the 
maximum recycling rate, Xo is the half-saturation constant, and q specifies the shape of the 
sigmoidal recycling curve. At s = l, r — 1 , xo = l, q — 8, and c = 0.53 (as in 

, the system has alternative stable states: a low phosphorous 
xl = 0.537, and a high phosphorous eutrophic state, xh — 1.491, 
separated by an unstable equilibrium (a saddle), xg = 0.971. 

The standard technique for studying systems like this one, is linear stability analysis. 
The eigenvalue of the linearized system at xg is Xs = 1.032, so it is an unstable equilibrium. 
The eigenvalues corresponding to xl and xh are Xl = ~ 0.899 and \h = ~ 0.797, respectively, 
so both xl and % are stable equilibria. The more negative the eigenvalue, the faster the 
return to the equilibrium following a small perturbation; \l<\h , so the linear analysis 
indicates that the oligotrophic state is more stable than the eutrophic state. 


Jayaprakash 2007 


oligotrophic state, 


Guttal and 


Ball-in-cup 

An alternative approach to quantifying stability, and one that is fundamental to the 
theory of alternative stable states, is the “ball-in-cup” heuristic (Beisner et al. 2003). In 
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this framework, the state of the system is represented by the position of a ball rolling on a 
surface. The ball rolls downhill, but is also subject to continual, stochastically varying 
perturbations. In the absence of perturbations, the ball will roll to the bottom of a valley. 
Such locations correspond to stable equilibria of the deterministic skeleton of the system 
(xl and xh in our example); a system with alternative stable states has more than one 
valley. The “cup” is the area surrounding an equilibrium that is attracted to it; this is 
called its domain (or basin) of attraction. 

The ball-in-cup framework is not just a useful metaphor - it can also yield a 
mathematical description. For the lake system, define 

U(x) = - [ X /(O^, (2) 

Jx H 


(£ is a dummy variable for integration), so that the differential equation becomes: 

^ = —U'(x). The dynamics of this system turn out to be equivalent to a ball-in-cup 
system with surface specified by the function U. In analogy with the physics of the ball-in¬ 
cup metaphor, U is called the “potential function” or simply the “potential”. For the lake 
system, this surface has local minima at xl and xh, as shown in figure [T)i. 

When random perturbations are present, the ball can be jostled from one basin of 
attraction to another. Note that stochasticity lies at the heart of the theory of alternative 
stable states. In a purely deterministic system, the ball would roll to an equilibrium and 
stay there. The presence or absence of other stable states would be irrelevant, because the 
ball would have no way of visiting them. Perhaps the surface could change over time, so 
that the basin of attraction occupied by the ball ceases to be a basin, and the ball rolls out 
to a different stable state. This situation corresponds to a bifurcation of the system’s 
deterministic skeleton; the ball’s transition requires the destruction of a stable state. In 
this paper, we are interested in how systems can transition between coexisting alternative 
stable states. Perturbations are required for the system to undergo these transitions; 
therefore, we argue that the appropriate framework for an alternative stable state model is 
a stochastic one. Furthermore, real ecological systems are always subject to random 
perturbations. In order to apply the ball-in-cup heuristic to a perturbed system, we next 
demonstrate an approach to incorporating stochasticity into model ([Tj). 


Stochastic Differential Equation Model 

If the nutrient concentration varies randomly over time, the lake can shift from one 
stable state to the other. To study this scenario, we translate the original deterministic 
model into a stochastic differential equation. A brief explanation of stochastic differential 
equation models is provided in appendix [Aj and more extensive accounts can be found in 
textbooks (e.g. Allen 2007). Here, we give an informal description of the major concepts, 
and use discrete-time analogies to avoid overly technical mathematical terminology. 

To emphasize that nutrient concentration is now a stochastic process, and not just a 
deterministic function of time, we switch notation from x(t) to Xft). For each i>0, x(t) is 
a number, but X(t) is a random variable, which can take on any of a set of possible values 
according to probabilistic rules. A realization of the stochastic process is a deterministic 
function of time associated with a specific set of random events; this can be thought of as 
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an observed time series, or the result of a single simulation run. 

In the original model Q, the external input of nutrients occurs at a constant rate c. In 
a small time interval dt, the external input is cdt. In reality, this input is likely to vary 
randomly; this is commonly modeled by adding a Gaussian white noise process, dW(t ) 
(“noise” is used synonymously with “stochastic” or “random”). At each £>0, dW(t) is a 
normally distributed random variable with mean zero and variance dt. Since the values are 
independent of t, this is simply written as dW. The white noise process we describe here 
has no temporal autocorrelation, and its frequency spectrum is uniform - the descriptor 
“white” is used in analogy with white light. The accumulated change obtained by adding 
dW over time yields a Wiener process, also known as Brownian motion. White noise is a 
useful starting point, but many applications require other types of noise; for example, 
colored noise might be used instead when perturbations are autocorrelated (e.g 


Sharma 


et al. 2014). A discussion about generalizing the framework in this paper to different noise 


types is included in the Limitations and Generalizations section. 

If the constant input rate c is perturbed by a Gaussian white noise process with 
intensity a (analogous to the standard deviation in discrete time systems), then the 
external input in a small interval dt is cdt + crdW. The change in nutrient concentration 
over this time interval is given by 


dX=\c-X + 


X q 


1 + x q 


dt + a dW. 


( 3 ) 


Again using equation (|2]) to define the potential, this system can equivalently be written as 

dX = -U\X)dt + a dW. (4) 


In terms of the ball-in-cup heuristic, the shape of the surface is specified by the potential 
function U, and this is independent of cr. The noise intensity cr only contributes to the 
movement of the ball on this surface, as determined by the last term in equation Q. 

We have described this model in terms of change over discrete time intervals, but it is 
also valid in the continuous time limit, dt— s-0. For continuous time, which will be the focus 
of the rest of this paper, ([3]) is called a stochastic differential equation. The notation in the 
stochastic differential equation dX =... is different than the deterministic differential 
equation notation ^ = ..., because the former must be defined using integral equations 
(the realizations of W(t) are not differentiable anywhere, so and hence would not 
make sense. We use the Ito integration scheme to define stochastic differential equations in 
this paper; see appendix |A| . 

Utility of the potential for understanding the stochastic lake eutrophication model 

One approach to understanding the stochastic lake eutrophication model is to calculate 
realizations (i.e. simulations) of (J3]) for particular values of cr. This approach is limited, 
because it requires setting a particular cr; we will see later that the potential function 
provides a more general way of studying system dynamics. A realization with a = 0.2 is 
shown in figure [1 ]d. All simulations in this paper were done with Mathematical and the 
code is available as a supplementary hie. The realization in figure [TJd, which is typical of 
realizations for this system with cr = 0.2, switches between the two stable states. It spends 
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more time near % than xl; this suggests that the eutrophic (higher phosphorous) state is 
more stable than the oligotrophic (lower phosphorous) state for this set of parameter 
values. Note that this behavior is in contrast to the results of the linear stability analysis of 
the deterministic skeleton. It is, however, in agreement with what the potential function 
tells us about the system, as we will demonstrate below. 

For ([3]), we find that U{xl) = 0.011, U{xs) = 0.047, and U(xh) = 0. Note that it is the 
relative, not the absolute, values of the potential function that are important, so the 
minimum value of the potential can be set at 0. U{xh) <U(xl), so the potential function 
indicates that the eutrophic state is more stable than the oligotrophic state. This 
corresponds to the intuitive notion that we obtained from examining realizations like the 
one in figure [TJd, but it contradicts the results from the linear stability analysis. This 
discrepancy arises because the linear stability analysis considers only an infinitesimal 
neighborhood of an equilibrium. In the presence of continuous stochastic perturbations, the 
system will leave such an infinitesimal neighborhood, and the linear analysis of the skeleton 
breaks down. The linear analysis provides information about the curvature of the potential 
surface at the bottom of basins of attraction, but this information is purely local, in that it 
does not take into account the larger geometry of the surface. Therefore, the potential 
function provides a more appropriate measure of stability for analyzing alternative stable 
states than linear stability analysis. 

The potential function also relates to other important features of the stochastic system. 
The probability density function, p(x,t), associated with the random variable X in (J3]) 
describes the probability that X{t) — x. It is the solution to the Fokker-Planck equation: 


dp(x,t ) 
dt 


d_ 

dx 


(U'(x)p(x,t)) + 


a 2 d 2 p(x,t ) 
2 dx 2 


The steady-state solution, p s (x) — lim p(x,t), is given by: 

t—> OO 


( 5 ) 



(-W )■ 


( 6 ) 


where Z — J 0 °°exp ( — 2 U a 2 ) dx is a normalization constant. This equation shows that the 


steady-state probability density is maximized at the values of x that minimize U , 
confirming that the minima (valleys) in U correspond to the most likely system states. 

The potential can be used to gain insight about the time it takes the system to switch 
between alternative stable states. If r5 s is the expected time it takes a trajectory starting 


at Xl to reach xh, (he., the mean first passage time), then (Kramers 1940): 


t Xh = 

X L 


7 exp I — (U(xs) — U (xl)) ) (1 + 0(a)). 

VU"(x L )\U"(x s )\ 1 M 


( 7 ) 


Swapping xh for xl yields a comparable expression for the expected time to reach xl from 
xh- The asymptotic notation O(-) describes the error of the approximation as a —> 0. The 
expected time for a trajectory to leave a basin of attraction around one of the stable states 
is thus largely dependent on the depth of that basin - the difference between peak U 
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(which occurs at the saddle equilibrium, xs) and the value of U at the stable equilibrium. 

The eigenvalue obtained in linear stability analysis describes the curvature of the 
potential at an equilibrium, equal to the second derivative of U\ it determines the prefactor 
that multiplies the exponential function in equation ([7]). For a fixed valley depth, increased 
curvature is associated with decreased mean first passage time. For instance, note that 
A l = —U"(xl)- As xl becomes more stable in the deterministic sense (i.e., as A l becomes 
more negative), the curvature at xl increases, and the mean first passage time decreases 
(similar statements hold for %). At first glance, this seems counterintuitive - increasing 
stability is associated with decreased escape time - but it makes sense because, for a fixed 
valley depth, increased curvature decreases the horizontal distance between equilibria. 

Knowledge about the potential function thus provides information about the steady- 
state probability distribution, mean first passage times, and transition frequencies, 
motivating its use as a stability metric (Zhou et al. 2012 Wang et al. 2011[ ). The potential 
function is especially useful because it does not depend on the noise intensity a (in contrast 
to the steady-state probability distribution and mean first passage times; see appendix [P|. 

Example 2: Consumer and Resource With Alternative Stable States 

If the potential is so good at quantifying biologically-relevant model behaviors, why isn’t 
it routinely applied in ecology? Unfortunately, in most cases, there will not exist a function 
U that satisfies the mathematical definition of a potential (see appendix [B]) . Systems that 
have such a function are called “gradient systems”. One-dimensional systems are always 
gradient systems, but systems with more than a single state variable almost never are. For 
non-gradient systems, we cannot use a potential function to quantify stability, as we did in 
the first example. It is for this reason that ecologists typically rely on approaches like 
linear stability analysis instead; although these approaches give more limited biological 
insights, they are more widely applicable mathematically. In what follows, we show how to 
generalize the potential for non-gradient systems, thus allowing us to apply the many 
desirable features of potential analysis to a much broader range of ecological systems. 

For an ecological example of a two-dimensional non-gradient system, we turn to a 
model of phytoplankton and zooplankton populations. Let R be the phytoplankton 
(resource) population density and C the zooplankton (consumer) population density. Using 


the deterministic skeleton of a standard plankton consumer-resource model (Collie and 


Spencer 1994, Steele and Henderson 1981), we obtain the stochastic differential equations 


dR = ( aR ( 1 ——^ I dt + o\ dW\ 

p J hi + R 2 


dC = 


pR 2 C 
k + R 2 


( 8 ) 


— pC 2 ) dt + (j-2dW‘2- 


Here W\ and W 2 are independent Wiener processes. The resource has logistic growth in the 
absence of consumers, with maximum growth rate a and carrying capacity /3. Consumption 
of resources is represented by a sigmoidal Type III functional response. 5 is the maximum 
consumption rate, and n controls how quickly the consumption rate saturates. 7 determines 
the conversion from resources to consumers. The consumers have a quadratic mortality 
term with coefficient p, which represents the negative impacts of intraspecific competition. 























<7 1 and cr 2 are the noise intensities for the resource and consumer populations, respectively. 

The additive form of the stochastic terms in this model represent random inputs and 
losses of resources and consumers. In situations where inherent growth parameters (e.g., a 
or 7 ) are stochastic, other forms of stochasticity would be appropriate. We will deal with 
additive noise here; the more general case is considered in appendix |F| 

We will analyze Q with parameters set at a = 1.54, /3 = 10.14, 7 = 0.476, 5 = k= 1, and 
/i = 0.112509. A phase plot of the deterministic skeleton is shown in figure [2^. The 
deterministic skeleton of this system has five equilibria: e 0 = (0,0), e^ = (1.405, 2.808), 
e B = (4.904,4.062), e s = (4.201,4.004), e P = (/3,0). 

A linear stability analysis shows that e 0 is an unstable equilibrium and e P is a saddle 
point, and e p are stable equilibria, and eg is a saddle point that lies between them. 
Equilibria and their stability are summarized in figure [2^. 

The eigenvalues of the Jacobian are —0.047 ± 0.4581 at and —0.377 and —0.093 at 
e p . For eA the real part of the eigenvalue with largest real part is —0.047, and for e p it is 
—0.093; therefore, the stability analysis concludes that e p is more stable, because this 
value is more negative than it is for 04 . 

A realization of the stochastic system (07 = 07 = 0.05, Figure^) shows switching 
between the two stable states. It is typical of most realizations we generated, in that it 
spends more time near (dotted white lines) than e p (dashed black lines). This 
realization, which had initial condition (xo, yo ) = (1, 2), spent 87% of its time in the basin 
of attraction corresponding to e^. Intuitively, it seems that should be classified as more 
stable than e p , but as in Example 1, this is not what was obtained via the standard linear 
stability analysis. 

Recall that realizations are of limited utility for stability analysis, because each value of 
a will produce different dynamics and different steady-state probability distributions (see 
appendix D and supplementary figure 1). The potential is defined independently of cr, and 
hence would be ideal for providing more general insights than cr-speci£c realizations. Of 
course, we do not have a potential function U for this or any other non-gradient system 
and hence cannot compare U(eA ) and f/(e p ). Instead, we turn to the Freidlin-Wentzell 
quasi-potential, which generalizes the notion of a potential. 

Generalizing The Potential 

For higher-dimensional models, we need to introduce a little bit of new notation. We 
can write an n -dimensional system of stochastic differential equations with additive noise 
as 

dX. = /(X) dt + a dW. (9) 

X= (Xi,..., X n ) is a column vector of state variables and W = (Ihj,..., W n ) is a column 
vector of n independent Wiener processes. We use the lowercase notation x= ( 27 ,..., x n ) 
to indicate a point in phase space (as opposed to a stochastic process). / is the 
deterministic skeleton of the system. It is a vector held: for every point x, /(x) specifies 
the direction that a deterministic trajectory will move, a is the noise intensity. More 
general ways of incorporating noise are considered in appendix [F} 

Following the same general approach as in example 1, the Fokker-Planck equation for a 


9 


( 10 ) 


two dimensional version of (J9|, with X = (X 1? X 2 ), x= (aq, x 2 ) and / = (/i,/ 2 ), is: 


dp 

dt 


d 
dx i 


(fiP) 


d 

dx 2 


(hP ) + 


2 


d 2 p d 2 p \ 

dx 2 dx 2 ) ' 


In the gradient case in Example 1, the steady-state solution of the Fokker-Planck equation 
was of the form (J6]) (replacing x with x and obtaining Z via integration over the positive 
quadrant). Here, there is no function U to play that role, but using the same general 
approach, assume that there is a function V (x) such that: 


p s (x) x exp ^- 2 ^ 2 X ^ • (11) 

The symbol x denotes logarithmic equivalence, details about which are in appendix |Bj 
When noise intensity is small, we can obtain an approximation for V (using asymptotic 
expansion; see appendix [D]) . This approximation, denoted by Vo(x), satisfies 

VVq ■ Wo + / • VF 0 = 0, (12) 


where the gradient operator V takes a scalar function ip as an input, and returns a vector, 
Vip = Jy), • • • j §,~~^ji that is the multi-dimensional analogue of the derivative. 

Intuitively, if one thinks of ?/;(x) as specifying the height of a landscape at a particular 
point x, then — V'0(x) points in direction of the steepest descent (as water would flow). 

Equation fll2| ) is the static Hamilton-Jacobi equation. Interestingly, Vo has key 
properties that make it a useful analog of a potential in a gradient system. First, Vo is 
independent of the noise intensity cr, just as the potential function U was in the gradient 
case. Second, if x(t) is trajectory of the deterministic skeleton of (J9]) , then 

j t (Vo (x(f))) = Wo ■ / (x(t)) = - W 0 ■ Wo < 0 , (13) 

and 4 (Vb (x(t))) = 0 only where VVq = 0. Thus Vo is a Lyapunov function for the 
deterministic system, which is an important feature for the ball-in-cup metaphor. If Vo(x) 
specifies an two-dimensional surface, then, in the absence of perturbations, trajectories will 
always move “downhill”. Again, this parallels the role that U played in the gradient 
systems. Third, we can interpret the relationship between / and the surface Vo. / is the 
deterministic skeleton that causes trajectories to move across the landscape, and - Wo is 
the component of / that causes trajectories to move downhill. The remaining component 
of /, which we denote by Q and call the “circulatory” component, is defined as: 


Q ( x ) = / (x) + VHq (x). 


(14) 


Vo satisfies the Hamilton-Jacobi equation, so Q ■ W 0 — f ■ W 0 + W 0 ■ Wo = 0, hence W 0 
and Q are perpendicular at every point. This motivates the label “circulatory” - in the 
absence of other forces, Q would cause trajectories to circulate around level sets of Vo- 
The function Vq generalizes the potential function to non-gradient systems and extends 
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to n-dimensional systems. Interestingly, Vo is a scalar multiple of a function called the 
Freidlin-Wentzell quasi-potential. The quasi-potential has extremely important properties, 
which we explore in the next section before applying all of these ideas to example 2. 

The Freidlin-Wentzell Quasi-potential 


Freidlin and Wentzcll (2012) analyzed stochastic differential equations using a large 


deviation principle, which is an asymptotic law determining the probabilities of different 
trajectories. These concepts can be best interpreted by imagining the state of the system 
(the position of the ball, or the current combination of population densities) being 
randomly perturbed within a “force held” imposed by the deterministic skeleton. Suppose 
the system starts at the stable state and travels to another state x. To complete this 
journey, the populations will need to do some “work” against the force held (i.e., they need 
to go “uphill”); this work is provided by random perturbations. Trajectories that require 
the least amount of work (require the least extreme stochastic perturbations) are the most 
likely. Suppose that 6{t ) specihes a path, parameterized by t, that goes from the stable 
equilibrium 0(0) = to another state 0{T) = x. T is total time it takes the populations to 
move along this path from to x. The amount of work required for the populations to 
follow a given path can be quantihed by a functional St called the action (see appendix [B] 
for details). 

In order to determine the amount of work it takes to get to some state x, one must 
minimize the action over all possible paths from c,\ to x, and all path durations T>0. The 
minimum action is called the quasi-potential, denoted <F eA (x). The quasi-potential depends 
on the starting point e^; when there are multiple stable states, the corresponding 


quasi-potentials can be stitched together to obtain a global quasi-potential, $(x) (Roy and 
Nauman 1995); see further details in appendix [Cj <3> is related to Vo by <h = 2 Vo 


(appendix [E]). In this paper, we use Vo instead of <3>, because Vo agrees with the true 
potential in gradient systems. The multiple of 2 in the relationship $ = 2 Vo is an 
inconvenient result of the Freidlin-Wentzell definition. Conceptually, these two functions 
measure the same properties, and computing one immediately yields the other. 

The quasi-potential can be calculated by solving the static Hamilton-Jacobi 
equation (12). This is a numerically difficult task, however; standard finite difference and 


finite element methods typically break down when applied to this kind of non-linear partial 


differential equation. Ordered upwind methods (Sethian and Vladimirsky 2001) are an 


innovative approach that circumvent the problems encountered by traditional methods. 

The basic idea is to create an expanding front of points where the solution is known, and 
march outward by considering and accepting solution values at adjacent points in ascending 
order. For use in systems of the form (|9]), the standard ordered upwind method was 
enhanced by Cameron (2012). Cameron’s algorithm allows for efficient computation of the 


quasi-potential. It forms the basis for QPot, a freely-available R package we have developed 


(Moore et al. Submitted) that includes a full set of tools for analyzing two-dimensional 
autonomous stochastic differential equations. QPot can be downloaded at CRAN [] or 
GitHub [^] To calculate the quasi-potential, users simply input the deterministic skeleton of 
the system, the domain, and the mesh size (although many other options are available). 


lr The Comprehensive R Archive Network, https://cran.r-project.org 
2 https: / / github.com/bmarkslash7 / QPot 
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Computation time for the ordered upwind method depends on the model and mesh size; 
example 2 took less than ten minutes on a fairly average personal computer. 

The Freidlin-Wentzell construction of the quasi-potential provides a mathematically 
rigorous justification for the Wentzel-Kramers-Brillouin (WKB) ansatz, which can be used 


to approximate mean first passage times in the small noise limit (Bressloff and Newby 


2014). The WKB method has been applied to calculate expected extinction times for 


several specific models in population dynamics and epidemiology (Meerson and Sasorov 

2009, 

Roozen 

1989, 

van Herwaarden and Grasman 

1995, 

Ovaskainen and Meerson 

2010) 


Example 2 Continued 

We generated solutions to the static Hamilton-Jacobi equation for the system Q using 
base points e^ and e B , and then matched them into a global quasi-potential by enforcing 
continuity at e^ and setting the minimum to 0. We divided this function by two to obtain 
Vq- The ordered upwind method was implemented using Cameron’s algorithm (Cameron 


2012). Mathematica was used for data processing and graphics generation, and the code is 


available as a supplementary hie. 

For the consumer-resource system ([8]), the resulting surface for Vo and a corresponding 
contour plot are shown in figure [3^-b. We find that Vo(eyi) = 0, Vo(es) = 0.007, 

Vo(ee) = 0.006. The relative values of Vo can be used to make calculations regarding first 
passage times and calculate transition rates between and e B . The most fundamental 
observation, however, is that Vo(e,4) < Vo(ee), which indicates that is more stable than 
e B . This contrasts with the linear stability analysis, but agrees with the qualitative picture 
obtained from realizations of the system. As in example 1, analyzing the system through 
the lens of a potential (or quasi-potential) function yields a completely different conclusion 
than the deterministic analysis, and one that aligns much more clearly with the simulated 
dynamics we observe. Furthermore, Vo(e,s) and Vo(e B ) are closer to each other than they 
are to Vo(e^). This indicates that e s and e B have similar stabilities, and it encourages us 
to move beyond the dichotomous classification of equilibria as either stable or unstable, 
which is often applied in linear stability analysis. The stable vs. unstable dichotomy 
classifies and e B as alike, and e B as different. The quasi-potential shows that it is e B 
and es that are alike, and that is different. By quantifying stability on a useful 
continuum, the quasi-potential offers a more nuanced perspective. 

Vo also provides a useful way to decompose the deterministic skeleton of equations ([8]) 
into physically interpretable parts, / = — VVo + Q. This decomposition is shown in 
figure [4ji-b. —VVo represents the part of the system that moves the system towards stable 
states, while Q represents the part that causes consumer-resource cycling. 

Example 3: Predator and Prey With A Limit Cycle 

The quasi-potential allows for stability analysis of attractors that are more complicated 


than equilibrium points. As discussed in Cameron (2012) and Freidlin and Wentzell (2012) 


and explained in appendix [Bj the quasi-potential can be defined for compact sets, such as 
limit cycles. As an example of a non-gradient system with a limit cycle, consider a 
stochastic version of the Rosenzweig-MacArthur predator-prey model (e.g. Logan and 































Wolesensky 2009): 


dR = [ aR [ 1 —— )-—— I dt + a\dW\ 

(3 J K~\~R 


( IRC 


dC = -— — /dC ) dt + a 2 dW 2 . 

\ k + R 


( 15 ) 


Here R is the resource density, C is the consumer density, and W\ and W 2 are independent 
Wiener processes. Consumption of resources is represented by a Type II functional 
response; otherwise the resource dynamics are the same as in example 2. In the absence of 
resources, the consumer density decreases at an exponential rate determined by p. oq and 
<j 2 are the noise intensity for the resource and consumer densities, respectively. We present 
the analysis of this model with a = 1.5, j3 = 45, 7 = 5, <5 = 10, k = 18, and fi = 4. 

Figure |2)o,d shows a stream plot of the system’s deterministic skeleton, and a realization 
with noise intensities oq = <72 = 0.8 over time interval [0,50]. This choice of noise intensity 
and time scale was made to illustrate clear population cycles with amplitude shifts. 

Surface and contour plots of V 0 for system (15) are shown in figure [3]>d. Recall that V 0 
provides a decomposition of the deterministic system into a “downhill” force and a 
“circulatory” force, as shown in figure [4]>d. In this case, — VVo causes trajectories to be 
attracted to the limit cycle’s trough. The circulatory component causes trajectories to 


cycle in this trough. This decomposition harkens back to Holling (1973), who made the 
following observation about dynamical systems: “There are two components that are 
important: one that concerns the cyclic behavior and its frequency and amplitude, and one 
that concerns the configuration of forces caused by the positive and negative feedback 
relations.” The latter is described by the gradient of Vo, the former by the circulatory 
component. Therefore, we see that the Freidlin-Wentzell approach provides a systematic 
way to distinguish between the two concepts identified by Holling. 

In this example, we cannot contrast the quasi-potential results with the traditional 
linear stability analysis, because the latter only applies to equilibrium points. 

Limitations and Generalizations 

In this paper, we have focused on applying the quasi-potential framework to stochastic 
differential equations models that share several characteristics: 1 ) time is continuous, 2 ) 
state variables are continuous, 3) noise is additive and the noise intensity is the same for 
both state variables, 4) noise is a direct perturbation to the state variables (as opposed to a 
perturbation to parameter values), 5) noise is white (as opposed to colored), and 6 ) noise 
occurs continually with low intensity (as opposed to occurring as discrete, abrupt events). 
For models with discrete state variables, different approaches in large deviation theory are 


needed (Wainrib 2013). However, our approach can be adapted to work in systems that 


deviate from several of the other characteristics. For instance, characteristic 1 is not a 
limitation of the quasi-potential framework; |Kifer (1990) describes how analogous concepts 
can be applied to discrete-time Markov chains (Kifer 1990 Faure and Schreiber 2014). 
Variable transformations (see appendix [F]) can be used to compute quasi-potentials for 
systems that deviate from characteristic 3 (e.g. those with noise terms of unequal intensity 
(oq 7 ^ cr 2 ), noise that scales with population density (demographic stochasticity; 
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(Ti\/JCidWi), or multiplicative environmental stochasticity ((TjAjdWj) (Hakoyama and 
Iwasa 2000)). Perturbations to parameters rather than state variables can be 


accommodated by explicitly modeling the parameter as a state variable with its own 
differential equation (Allen 2007). A similar approach can be applied to models with 
colored noise (i.e., models that do not have characteristic 5). The noise process itself can 
be explicitly modeled as a state variable with its own differential equation (e.g., an 
Ornstein-Uhlenbeck process). Unfortunately, increasing the dimensionality of the state 
space in these ways makes the process of numerically calculating the quasi-potential even 
more challenging. Given the pace of development of numerical techniques (Cameron 2012), 
however, it is conceivable that solving such systems will soon be more practical. 

Characteristic 6, which states that noise occurs continually with low intensity, is central 
to the quasi-potential framework. The expressions relating the quasi-potential to 
steady-state probability distributions and mean first passage times are based on the 
assumption that the noise intensity is very small. As a rule of thumb, these approximations 
are only useful when a 2 is much less than 2 AVo, where AVo is the difference in the 
quasi-potential between the stable equilibrium and the saddle. In appendix [H] we provide 
details on how mean first passage time scales with noise intensity, and present a numerical 
examination of these concepts applied to example 2. For systems that experience extreme 
events and external shocks (e.g., natural disasters, extreme climactic conditions, invasive 
species introductions, etc.), the quasi-potential no longer provides complete information. If 
a shock directly impacts the state variable (e.g., if the lake system in example 1 were to 
receive a massive pulse of phosphorous run-off), the ball in the ball-in-cup diagram would 
experience a large, instantaneous horizontal displacement (perhaps skipping over 
intervening valleys and hills). If the system reverts to deterministic dynamics, or stochastic 
dynamics with lower-intensity perturbations after the shock, the quasi-potential will still be 
useful for describing the system’s response after the shock. In the presence of large shocks, 
though, the quasi-potential loses its ability to make probabilistic predictions. If a shock 
impacts the state variable indirectly (e.g., if an invasive species entered the lake and 
fundamentally altered the phosphorous cycling), the shape of the quasi-potential surface 
would change dramatically. The interaction between a dynamically changing quasi¬ 
potential surface and state-variable noise would be difficult to analyze using the methods 
presented here. 

The three examples in this manuscript show that the quasi-potential often provides a 
more informative stability metric than traditional linear analysis. Linear stability is much 
easier to measure in the held, though. This can be done by slightly perturbing a system 
and measuring the time it takes to return to equilibrium. Before the quasi-potential can be 
calculated, a model must be fit to observed data and validated. This limitation is also 
shared by other methods for analyzing systems with alternative stable states, which 
depend explicitly (e.g. Boettiger and Hastings 2012) or implicitly (e.g. Dakos et al. 2008) 
on underlying models. Fortunately, carefully controlled experiments (Dai et al. 2012) and 
advances in model-fitting (Ives et al. 2008) point toward a promising future for the 
empirical study of shifts between alternative stable states through models. 
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A Path Through the Quagmire of Stability Concepts 

Systems with alternative stable states are only interesting when perturbations can 
cause shifts between states; when these stochastic perturbations are continual and random, 
as in most ecological systems, stochastic models are appropriate. When state and time 
variables are continuous, stochastic differential equations like ([9]) are the best option. The 
three examples presented in this paper show that the quasi-potential provides a useful way 
to study such stochastic differential equation models. In particular, it provides a way to 
quantify the relative stability of alternative stable states. 

Unfortunately, many notions of stability were developed for a deterministic context, 
and these can be misleading when applied to stochastic systems (as in examples 1 and 2). 


Our goal is not to add to the existing tangle of stability definitions (Grimm and Wissel 


2008), but rather to provide a clarifying mathematical interpretation. Many existing 


definitions can be related to the ball-in-cup heuristic, and the quasi-potential shows that 
this metaphor has a useful and rigorous mathematical meaning. The translation between 
mathematical model and potential surface is easy in gradient systems (in particular, for 
one-dimensional systems, which are always gradient systems). The translation for more 
general systems is less obvious, but the quasi-potential fills that need. 

Figure [5^ is a ball-in-cup diagram of the potential for a one-dimensional system that 
helps to illustrate several important concepts associated with stability. These concepts are 
equally relevant for higher dimensional systems, where the ball rolls on a multi-dimensional 
surface specified by Vo (half the Freidlin-Wentzell quasi-potential) instead of a curve. 

One metric of stability for an equilibrium e 0 is the curvature of Vo at e 0 (dashed black 
line in figure^). The greater the curvature, the more difficult it is to perturb the system 
away from eo, and in this sense, the more stable eo is. In one dimension, the curvature at 
eo is U"(eo), which is minus the eigenvalue obtained in linear stability analysis. In higher 
dimensions, the eigenvalues are again directly related to curvature, now along different 
planar sections of Vo (see appendix [G]) . Thus, measuring the curvature of Vo at e 0 is 
equivalent to determining asymptotic stability through linear stability analysis. 


Asymptotic stability has a long history in ecology (May 1973). The primary problem 
with this metric is that it is purely local - once a trajectory is perturbed outside of a tiny 
neighborhood of an equilibrium, nonlinear effects can come into play and the 
approximation is no longer informative. Furthermore, this approach views perturbations as 
being isolated one-time events. With this view, a system is displaced, and then the 
dynamics proceed deterministically without further perturbation. In reality, perturbations 


often take place on a continual basis. Indeed, as noted by Ives (1995), “To apply generally 


to ecological communities, stability needs to be defined for stochastic systems in which 
environmental perturbations are continuous and equilibrium densities are never achieved.’ 


Likewise, Neubert and Caswell (1997) write, “real ecosystems are seldom if ever subject to 


single, temporally isolated perturbations. Nevertheless, our analyses, together with most 
theoretical and experimental studies of resilience, ignore the effects of continual stochastic 
disturbances in the hope that the deterministic results will shed light on the stochastic 


case. 


A second metric of stability of an equilibrium e 0 is the minimum distance between e 0 
and the boundary of its domain of attraction (dotted line in figure [5^). The width of the 
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basin of attraction measures the magnitude of perturbation that a system can sustain and 
still be guaranteed to return to eo- One problem with this metric is that, like asymptotic 
stability, it views perturbations as singular, isolated events. For this metric, it is only the 
boundary of basins of attraction that matter, not the shape or height of Vo- If 
perturbations happen continuously, the shape and height of the Vo are important. 
Nonetheless, this basin width metric can be extremely useful. 

A third metric of stability is the height of Vo (gray line in figure [5^,). Holling (1973) 
anticipated this concept, and called it resilience, which he explained with ball-in-cup 
diagrams. He defines one aspect of resilience, writing: “the height of the lowest point of the 
basin of attraction ... will be a measure of how much the forces have to be changed before 
all trajectories move to extinction of one or more of the state variables”. Holling had no 
way of defining the surface, and so could not actually quantify notions like “height”; the 
quasi-potential solves this problem. Holling’s identification of the difference between 
asymptotic stability and this definition of resilience (basin height) is hugely important, and 
it has major consequences for the analysis of alternative stable states. 

This third metric is perhaps the most useful of the three we have explored. Unlike the 
first two metrics, it is appropriate for use in systems that undergo continuous stochastic 
perturbations. As we saw in the examples in this paper, it can be used to compute mean 
first passage times, and is directly related to steady-state probability densities. 

These three metrics of stability can yield conflicting information about alternative 
stable states. Figure |5 Jd shows these three metrics for the equilibria and e# from 
example 2. Note that the basin width metric and the quasi-potential metric show that 
is more stable than e#, but the asymptotic stability metric shows the reverse. Appendix [T] 
demonstrates that the equilibria in a multi-stable system can exhibit any combination of 
the three stability metrics. That is, one equilibria can be classified as most stable according 
to the first metric, but not the second or third; or by the first and second, but not the 
third; etc. 

Resilience is a concept closely related to stability, and like stability, it is defined in 
different ways by different authors. In a large review of the ecological literature, 


Myers-Smith et al. (2012) found that resilience was used in many ambiguous and 


contradictory ways. Some authors, like Holling (1973) view stability and resilience as 
distinct properties; others, like Harrison (1979) define resilience as a single aspect of 
stability. Pimm (1984) and|Neubert and Caswell (1997) define resilience as essentially the 


asymptotic stability metric, while Harrison (1979), Peterson et al. (1998), and Gunderson 


(2000) define it as essentially the basin width metric. Ives and Carpenter (2007) defines 
Holling’s resilience using the dominant eigenvalue of the saddle that separates alternative 
stable states; like the asymptotic stability metric, this is the result of applying a local 
analysis to the deterministic skeleton of a system. 


Hodgson et al. (2015) argue that resilience cannot be quantified by a single metric, and 


use a potential function to illustrate the different components of resilience, which include 
latitude (the width of the basin of attraction) and elasticity (the asymptotic stability 
metric). The quasi-potential framework aids this clarification about resilience by extending 
it to multi-dimensional systems. 

The quasi-potential is also useful for understanding several other concepts related to 
stability. Reactivity (Neubert and Caswell 1997) differs from asymptotic stability, in that it 
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quantifies the immediate (as opposed to long-term) growth or decay of perturbations. In 
the quasi-potential framework, reactivity is related to the circulatory component of the 
vector held. In the neighborhood of asymptotically stable equilibria with high reactivity, 
the circulatory component of the vector held will carry trajectories away from the 
equilibrium before bringing them back. 


Harrison (1979) defined resistance as the ability of a system to avoid displacement 
The stress is quantified in terms of an environmental parameter 


during a time of stress, 
distinct from the state variables, and hence the interpretation of resistance depends on the 
parameter under examination. Resistance is best viewed as a measure of how dramatically 
Vo changes due to environmental parameter changes. 

Finally, Harrison dehned persistence as the ability of a system to stay in a given range 
when continual perturbations are applied. He notes that this is the property that is most 
biologically useful, and that stochastic differential equations are the best mathematical 
modeling tool to assess it. Unlike his definitions of resilience and resistance, this definition 
views the dynamics of the system as stochastic and subject to continual perturbations. He 
was unable to venture far with the mathematical analysis for this definition, but the quasi¬ 
potential provides a way forward. Mathematically, persistence can be dehned as the first 
passage time for a system to leave a specified domain, which is directly related to the 
quasi-potential. Thus Harrison’s persistence is another manifestation of the quasi-potential. 

Despite the confusing array of stability concepts currently used in ecology, we believe 
that the quasi-potential concept provides hope for clarity. The three metrics associated 
with the quasi-potential show how many of these concepts are deeply related (figure [5]). 

The mathematics developed by Freidlin and Wentzell (2012), coupled with numerical 
advances by Cameron (2012), make the quasi-potential a practical and accessible tool for 
ecologists to study alternative stable states. This paper’s goal is to demonstrate the utility 
of the quasi-potential, and to properly position it in terms of existing ecological ideas. 
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Figure 1 . Lake eutrophication model (example 1). (a) The potential function for equation ([Tj). 
The horizontal axis is the scaled nutrient (phosphorous) concentration and the vertical axis is 
the (dimensionless) potential. Gray disks are stable equilibria, and the white disk is an unstable 
(saddle) equilibrium. The dynamics of the system can be represented as a ball rolling on the surface 
specified by the potential function. Note that the basin around xh is deeper than that around xl- 
(b) A realization of equation ([3]), which models nutrient concentration, x, as a function of time, 
t. Variables are scaled, so the units are dimensionless. Integration was performed with the Euler- 
Maruyama method and At = 0.005. The solid lines corresponds to stable equilibria xl (lower) and 
xh (higher) for the deterministic skeleton. The dashed line corresponds to the the saddle point xs 
of the deterministic skeleton. Note that the realization spends more time near xh than near xl- 




21 

























Figure 2. (a) Stream plot for the deterministic skeleton of the consumer-resource model in 

example 2. Unstable equilibria are white disks and stable equilibria are gray disks. The unstable 
equilibrium ep is not shown, but would appear on the x-axis to the right of where the graph 
is truncated. Variables are scaled, so the units are dimensionless. Lines and arrows show the 
direction of trajectories for equations ([8]) in the absence of noise. (b) Similar stream plot for 
the deterministic skeleton of the consumer-resource model in example 3. The white disk is an 
unstable equilibrium and the gray line is a stable limit cycle, (c) A realization of equations Q 
for example 2. Integration was performed with the Euler-Maruyama method and At = 0.025. 
Resource population density is black and consumer population density is gray. The dotted white 
lines correspond to the equilibrium e^, and the dashed black lines to the equilibrium ep. (d) A 
realization of equations (15) for example 3, with a = 0.8. Integration was performed with the 
Euler-Maruyama method and At = 5 x 10 -4 . Resource population density is black and consumer 
population density is gray. 
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Figure 3. (a) The quasi-potential function for the consumer-resource model, equations ([8]). Vari¬ 
ables are scaled, so the units are dimensionless. Note that the quasi-potential surface is much deeper 
around than es- The quasi-potential is truncated at 0.02 for display purposes; it continues to 
increase in the regions outside the plot, (b) Contour plot for the same model. The white disk is 
the saddle point eg. The gray disks are the stable equilibria and e#. (c) The quasi-potential 
function for equations (15). (d) Contour plot for the same model. The white disk is an unstable 
equilibrium, and the white dashed line is a stable limit cycle. 
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(a) (b) 
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Figure 4. (a) and (b) are the orthogonal decomposition of the deterministic skeleton of the 

system Q. (a) The “downhill” component, — VVo- (b) The “circulatory” component, Q. Gray 
disks are stable equilibria. The white disk is an unstable equilibrium, (c) and (d) are the orthogonal 
decomposition of the deterministic skeleton of the system ([8]). The thick gray line is a stable limit 
cycle. 
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Figure 5. (a) A schematic diagram of the relationship between various concepts of stability, as 

related to the quasi-potential and Vo- (b) A comparison of three different metrics of stability for 
the system ([8]) . 
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Appendices 

A Stochastic Differential Equations 

In example 1, we briefly described a stochastic differential equation model for lake 
eutrophication. In this section of the appendix, we provide background information about 
stochastic differential equations. A random variable X is a variable whose value is subject 
to chance. When a specific outcome X = x is observed, it is called a realization. A 
stochastic process X(t) is a family of random variables indexed by the parameter t, which 
usually represents time. Time can be measured discretely or continuously; this latter case 
falls in the realm of stochastic differential equations. A realization, X{t) = x(t), is obtained 
when the stochastic process is observed at each time t. Note that a realization x{t) is a 
deterministic function of time. 

A continuous-time stochastic process of particular importance is the Wiener process, 
also known as Brownian motion, and denoted by W{t). This process can be visualized as 
the limit of a discrete time random walk, which changes by an amount AW per each time 
step At. Each increment AW is selected from a normal distribution with mean 0 and 
variance At. The Wiener process is the limit of this random walk as At —> 0. It turns out 
that the Wiener process is completely characterized by three properties: 

1. W(0) = 0 

2. W ( t ) is almost surely continuous everywhere. This means that, with 100% 
probability, a realization will be continuous (aside from possibly a few bad points, 
which have measure zero). 

3. If 0 < si < t\ < S 2 < t 2 , then W{ti) — W(s i) is normally distributed with mean zero 
and variance t\ — .si, W{t- 2 ) — W(s 2 ) is normally distributed with mean zero and 
variance t 2 — s 2 , and W(ti) — W(s i) and W(t 2 ) — W(s 2 ) are independent. 

The Reimann-Stieljes integrals of elementary calculus are defined as the limits of finite 
sums. Integration with respect to a Wiener process can be defined in a similar way. The Ito 
integral of the function h of a stochastic process X(t) over the interval [0, T] is defined as: 

pT n -1 

/ h(X(t)) dW = lim V h{X(ti)) ( W(t i+1 ) - W(U)), (16) 

I n n—>oo L —* 

J{J i =0 

where {[£*, tj+i)}” =0 is a partition of [0, T\. Note that this integral is a stochastic process 
itself; each realization of X and W leads to a different realization of the integral. In the Ito 
integral, h(X(t)) is evaluated at the left end points of the intervals of the partition. If a 
trapezoidal rule is used instead, then the result is the Stratonovich integral. In this paper, 
we use the Ito integral, because of the way it discriminates between the past and the 
future. A process X(t) is called “non-anticipating” if its value at t is independent of values 
of W(s), for s > t. If X(t) is non-anticipating, then the Ito integral defined above is, too. 
The Stratonovich integral is not, because calculating the integral at time s, ti < s < t* + i 
requires knowledge of X(ti + 1 ). Basically, the Ito integral cannot “see into the future”, 
while the Stratonovich integral can. 

Having defined integration with respect to a Wiener process, we can now define a 
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stochastic differential equation. Consider a deterministic autonomous differential equation, 

I = '<*>■ < 17 > 

In a small time period At, the variable x changes by an amount of approximately 
Ax = f(x) At. Now suppose that the variable x{t) is subject to random disturbances, and 
hence is a stochastic process X(t). To approximate the value of this stochastic process at 
time T, we discretize time into m small intervals, each of length At. Let Ah = X (iAt), and 
AXi = X l+ \ — X j. During a time period of length At, there are probably many small 
perturbations that affect A"; if they have finite variance, then by the central limit theorem, 
adding these small perturbations up yields a normally distributed random variable. We will 
assume that this accumulated perturbation over a time period of length At has mean 0 and 
variance a 2 At (the linear relationship with At is required in order for X(T) to have finite, 
non-zero variance in the continuous time limit). Therefore, the change in the stochastic 
process over a time interval of length At can be written as: 


AXi = /(Ah) At + aAWi, 


(18) 


where Ah) is normally distributed with mean 0 and variance At. Adding up the changes 
in the process over the time interval [0, T] yields 


X(T) = A’10) +'£,J(X i )At + aJ2 AW„ 


(19) 


2=1 


2=1 


which suggests an integral equation for the continuous time limit, 


X(T) = X(0) + / f(X)dt + a dW. 


( 20 ) 


If the intensity of perturbations depend on the value of X, then equation (20) can be 
generalized to 

X (T) = X(0) + [ f{X)dt + a [ g (A) dW 


( 21 ) 


The integrals in equation (21) make the notation cumbersome. In light of this, a modified 
notation is used. The stochastic differential equation 


dX = /(A) dt + a g(X) dW 


( 22 ) 


formally means that X(t) is a solution to equation (22). 


Note that does not exist, 


because the sample paths of W(t) are almost surely nowhere differentiable. This is why the 
notation in equation (22) is used; it reminds us that A (t) is defined by the integral 
equation (21). 
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B Freidlin-Wentzell quasi-potential 

In this section, we provide a more formal definition of the Freidlin-Wentzell quasi¬ 
potential. Consider a system of stochastic differential equations 


dX = /(X) dt + a g(X) dW, 


(23) 


where X = (Ad, ..., X n ) is a vector of state variables, W = ( W \, ..., W rn ) is a vector of 
m independent Wiener processes. Vectors in this paper should be interpreted as column 
vectors. The lower-case notation x = (aq, ..., x n ) is used to indicate a point (as opposed to 
a stochastic process). / is a vector field that is the deterministic skeleton of the system. 
g(x) is a matrix that determines how the different noise sources affect the state variables, 
and a is the noise intensity. For simplicity, we will focus on the case where m = n and g(x) 
is the identity matrix, which represents constant-intensity isotropic noise, affecting each 


state variable with equal intensity. Under these assumptions, equation (23) can be written 

(24) 


as 


dX = f(X)dt + adW. 

In appendix [F| we will return to the general case (23), but constant, isotropic noise 
provides a useful starting point. If there exists a function U(x) such that / = — W, then 
the differential equations are called a gradient system, and the function U is called a 
potential function. Like one-dimensional systems, a multi-dimensional gradient system can 
be viewed with the ball-in-cup framework. For n — 2, the relevant metaphor is a ball 
rolling on a two-dimensional surface specified by the function C/(x). For n > 3, the 
situation is difficult to visualize, but the same general intuitive aspects hold. The steady- 
state probability distribution of higher-dimensional gradient systems is related to the 
potential U in the same way as in ([6]), except x replaces x and Z is obtained from an 
n-dimensional integral. Expressions for the mean first passage time between stable 
equilibria separated by a saddle are similar to the one-dimensional case as well. 


Unfortunately, gradient systems are a very special situation. In most cases of (24), there 
will not exist a function U satisfying / = — W. For these non-gradient systems, we cannot 
use a potential function to quantify stability, as we did in example 1. In what follows, we 
develop an approach that is conceptually analogous but applicable to non-gradient systems. 

In the following, we will use the concept of logarithmic equivalence, denoted by x. We 
write f(x) x e Kh <*) if 

lim /W 1 In (f(x)) = h(x). (25) 

ft—>-00 

The Freidlin-Wentzell approach is to obtain a large deviation principle for trajectories x(t) 


of (24). In this context, a large deviation principle is an asymptotic rule that determines 


how likely it is for realizations of (24) to depart from a given path. To make this concrete, 
let a be an asymptotically stable equilibrium of f in (24). Let bel" and T > 0. Let ©t 
be the set of all absolutely continuous paths 6 : [0, T] —y R n such that 0(0) = a and 


6(T) = b. We will study the probability that a realization x CT (t) of (24) with noise intensity 


a and with x cr (0) = a and x CT (T) = b stays close to 6 e Q T . A large deviation principle 









declares that there exists a <5 0 > 0 such that, if 0 < 5 < £ 0 , then 


( 


Pr < sup |x, 

lo <s<T 


'(s)-6(s )| < 5 


exp 


S T (0) 


cr z 


(26) 


where the logarithmic equivalence holds as a —> 0. The functional St ■ ©t —t [0, oo) is 
called the action, and it is defined by 


St{9) = 2 


fm) - m 


dt. 


(27) 


Note that St measures how much 8 deviates from the vector held /. If St(&) = 0, then 8 is 
a trajectory of the deterministic system, ^ = f(x). The action St is related to the 
probability distribution of X by 

lim a 2 In (Pr (X(T) e fi|X(0) = a}) = - inf (S T (0)|0(O) = a ,8(T) e ft}, (28) 

cr —>0 9£@t 

where hi is a domain in M n . For details on the technical assumptions behind this 
relationship, see Freidlin and Wentzell (2012). To get from a to b in a “likely” way, the 


action should be made as small as possible. This motivates the definition of the Freidlin- 
Wentzell quasi-potential (or simply quasi-potential), $ a : M n —» [0, oo), 


*a(b) = mf n (5 r (0)|0(O) = a ,8(T) = b} . 

1 >U,C7Ev7X' 


(29) 


Note that the inhmum is taken over paths of all durations (that is, all times T > 0). 

The quasi-potential is the value of the action for the minimum-action path (i.e., the 
most likely path) between a and b. It is closely related to first passage times from domains 
of attraction. If D is a region contained within the domain of attraction of a, then the 
expected time until a trajectory exits D, t® d , is given by 


l™M t» d ) = inf *.(x). 

CT—>-0 x£oD 


(30) 


The quasi-potential need not be defined solely in terms of an isolated asymptotically 
stable equilibrium a. Cameron (2012) generalized the quasi-potential, and defined it for 


compact sets. This generalization allows the quasi-potential to be determined for limit 
cycles (as demonstrated in example 3). A different approach to generalizing the 


quasi-potential to compact sets can be can be found in Freidlin and Wentzell (2012). 


Cameron’s generalization requires considering the geometric action (Heymann and 


Vanden-Eijnden 2008&fa ), which we will denote by S*. Suppose that 8 G ©t, and ^{y) is a 
reparameterization of 8 such that ^(0) = 8(0) and V’(r'o) = 8(T). Then the geometric 
action is 


-v 0 


s*W0= / l/W>M)l WOI - fbK")) ■4’(")dv. 


( 31 ) 


The value of S* is independent of the parameterization of ijj. If A and B are compact sets 
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in M n , then the quasi-potential can be defined by 


$a(B) = inf {S’(4’)\4’(0) e A,4>(v o) 6 B} . 


(32) 


C A global quasi-potential 

In systems with multiple stable equilibria, it is desirable to obtain a global quasi¬ 
potential that describes how trajectories switch between states. In the preceding section, 
the quasi-potential was defined in terms of a stable equilibrium a. Suppose now that there 
are two stable equilibria, ax and a 2 , with corresponding domains of attraction Zifi and D 2 . 
The action functionals can be used to obtain <h ai and $ a2 , but these quasi-potentials are of 
limited utility outside of D\ and D 2 , respectively. The minimum action path from ai to a 2 
will follow streamlines of the vector field once it enters Do. This will result in no 
accumulated work; hence <F ai will be flat along streamlines in D 2 . The quasi-potentials 
both describe dynamics well within their domains of attraction, but in order to create a 
complete surface in the spirit of a classical potential function, it is necessary to combine 
the two. This is easily accomplished if there is a single saddle point s that lies on the 
separatrix between D\ and D 2 . We find the constant 

C = $ ai (s)-«h a2 (s) (33) 


so that = <h a2 + C agrees with <F ai at s. Finally, we compute the global quasi-potential 
<f> as 


4>(x) = min (* ai (x),4>;,(x)) . 


(34) 


More complicated cases can arise when domains of attraction are connected by more than 
one saddle (Freidlin and Wentzcll 2012). For details about how to combine local 


quasi-potentials into a global quasi-potential in these more complicated cases, see Freidlin 


and Wentzell 

(2012) 

Moore et al. 

(Submitted 

), and 

Roy and Nauman 

(1995) 


D Small noise expansion of V 


This section describes the relationship between V and Vo, and shows the derivation of 
the Hamilton-Jacobi equation for Vq. The Fokker-Planck equation associated with the two- 


dimensional version of ( 

24 

) is ( 

10 

). Under relatively mild conditions on the function / (for 

details, see 

Freidlin and Wentzcll 

2012 

), there will exist a steady-state probability 


distribution 


p s (x) = lim p(x,t). (35) 

t—> OO 


Steady-state distributions can often be approximated by very long-time realizations. 
Determining when such an approximation holds is the subject of ergodic theory (see 


Arnold 2010). Approximations to steady-state distributions for the consumer-resource in 
the system 


example 2 (i.e., the system (|8j) ) are shown in figure SI. Each panel corresponds to a 
different noise intensity, a. Qualitatively, figure SI confirms that trajectories spend more 
time near than e#, and hence a sensible stability metric should classify as more 
stable than e#. However, it also clearly shows that the steady-state distribution depends 
on the noise level; each choice of a yields a different distribution. If one is interested in the 
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general properties of the system, and not just the steady-state distribution for a specific 
noise intensity, then steady-states distributions are of limited utility. 

The “effective potential” (not to be confused with the potential or quasi-potential) is 
defined as 


l/(x) 



lnps(x) + C , 


(36) 


where C is a constant. The effective potential’s relationship with the steady-state 
distribution makes it a helpful tool. The peaks of the steady-state distribution correspond 
to valleys of the effective potential, and vice versa. There are two reasons why we do not 
adopt the effective potential as a stability metric in this paper. First, the effective potential 
depends on a, and hence suffers from the same issue as the steady-state distribution. In 
the ball-in-cup metaphor, the noise intensity a determines the perturbations of the ball as 
it rolls, rather than determining the shape of the landscape. Second, the effective potential 
is not a Lyapunov function for the deterministic system, so a trajectory of a system with 
zero noise does not necessarily move downhill. Finally, a decomposition based on the 
gradient of V is not orthogonal. Despite these shortcomings, the effective potential is 
closely related to the quasi-potential. Solving (36) for p s (x) yields 


2 C 2 V(x) 

p s (x) = e^e 


(37) 


Substituting this into the Fokker-Planck equation yields 


dV 


dV a 2 


iwr + f 1 — + f 2 —- — iv 2 v+ 


dx\ 


dxo 


9fi | dh 

dxi dx 2 


= 0. 


(38) 


To simplify this equation, we consider how the system behaves for small noise values, and 

2 

expand V in terms of the small parameter e = This yields 

OO 

r(x) = £ v,( 3 


x e 


(39) 


i=0 


where V) is the coefficient function associated with order e*. Inserting this into (38) and 
retaining lowest-order terms, we obtain the Hamilton-Jacobi equation for Vq 


|VDo| 2 + /i^ + / 2 ^=0. 


dx\ 


J dx 2 


(40) 


E Hamilton-Jacobi equation for the quasi-potential 

By deriving the Hamilton-Jacobi equation for the quasi-potential, we can verify the 
relationship <3> = 2Vo- This relationship is crucial. We described key properties about Vo 
concerning the effective potential, the steady-state probability distribution, the Lyapunov 
property, and the orthogonality of the decomposition / = — W + Q ; in appendix B, we 
described properties of $. The relationship $ = 2 Vo shows that these functions share those 
properties; they only differ by multiplication of a scalar. 

Bellman’s Principle from optimal control theory can be used to derive the Hamilton- 
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Jacobi equation for $. We sketch the proof from Cameron (2012). The calculation of 
$,4 (x), the value of the quasi-potential starting at a compact set A and going to a point x, 
can be viewed as an optimal control problem. We seek to minimize the value function 
$a(x) by choosing an optimal path This path is controlled by the velocity vector 

ip{v). We are free to choose the parameterization of ^(u), so we select one where the 
velocity vector has unit magnitude at every point. The optimal control problem amounts 
to determining the tangent direction ip{v) for each i/, so that the resulting path minimizes 
the action. Bellman’s Principle essentially turns this problem into a recursive equation. 
Heuristically, one can imagine the last segment of an optimal path ^{v) from ^(0) G A to 
= x. This last segment is specified by the parameter values v G [K — 5, K]. Clearly 
this optimal path will be optimal over the interval [K — 5, K]. Therefore, if one knows the 
optimal path up to parameter value K — 5 , one knows the remainder of the path as well. 
Mathematically, this principle takes the form: 


*K 


$a(x) = inf 

beS " -1 


A small 5 expansion yields: 


Iks 


\f{^{v))\ - /(V’O)) • V’W du + $A (ip(K - 5)) 


$a(x) = inf (|/(x)| -/(x) • ^)8 + $a(x) - V$a(x) • ip 8 

be S ’ 1 — 1 


(41) 


(42) 


inf 

be S 71-1 


Solving this equation is equivalent to solving: 

l/( x )l ~/( x ) •'0- V$ A (x) -^} =0. (43) 

Using the Cauchy-Schwarz inequality, one finds that the infimum of the left-hand side 

/(x) + V$a(x) 


of (43) occurs when 


= 


l/( x )l 


Substituting this into (43) yields 


|V$a| + 2 V$a • / = 0. 


(44) 


(45) 


Comparing this to (12), we can see that, if solutions exist, they have the relationship 


$a = 2 Vq- Classical solutions do not always exist for the Hamilton-Jacobi equation, so it is 


often necessary to consider a class of weak solutions called “viscosity solutions” (Sethian 


and Vladimirsky 2001, Crandall and Lions 1983, Crandall et al. 1984). When a classical 


solution does exist, it coincides with the viscosity solution. 
F Other noise structures 


This paper focuses on the case g(x) = I in equation (23), where / is the identity 


matrix. The quasi-potential can be calculated for more general cases. Such a generalization 
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requires a modification in the definition of the action: 


r-T 




(46) 


h3 


where q(x ) = ( g(x)g T (x )) 1 . The large deviation relationships, (26) and (28), are still 
valid. The Hamilton-Jacobi equation for (23) is 


E 




Qi. 


d<$> 5$ 

3 dxi dxi 


2 V<3> ■ / = 0. 


(47) 


Alternatively, one can find a transform of (23) that turns the system into the form (24), 


compute the quasi-potential in these new coordinates, and then back-transform to the 
original coordinates. For the system 

dX x = ,/i (X,, X 2 ) dt + cr gi dWi 
dX 2 = f 2 (X u X 2 ) dt + ag 2 dW 2 , 


(48) 


(49) 


where g± and g 2 are constants, the appropriate transform is Ad = g 1 1 Ad, Ad = g 2 1 Ad. For 
the system 

dX d = /i (Ad, Ad) dt + a g x Ad dW 1 
dX 2 = f 2 ( Ad, Ad) dt + a g 2 X 2 dW 2 , 

the appropriate transform is X\ = g 7 1 ln(Ad), X 2 = g^ 1 ln(Ad). 

G Curvature 

The concept of curvature is more nuanced for surfaces than it is for curves. The 
principal curvatures of the surface specified by Vo at e 0 are the largest and smallest 
curvatures of the one-dimensional normal sections at e 0 . A normal section is obtained by 
intersecting a plane containing the normal vector of the surface Vo at e 0 with Vo. The 
principal curvatures correspond to the eigenvalues of the Hessian matrix of Vo. In the 
gradient case, / = — VVo, so the Hessian matrix of Vo is simply the negative of the 
Jacobian matrix of /. I 11 other words, the principal curvatures of the surface Vo are the 
eigenvalues obtained from the linear stability analysis (except with the sign changed). 

H Mean first passage time asymptotics 

Steady-state probability densities and mean first passage times can be determined from 
the quasi-potential, but only in the small-noise limit. These quantities are often expressed 
in terms of logarithmic equivalence as a —> 0. Accurate calculation involves not just the 
exponential part of the relationship, but also the prefactor. For a gradient system with 
potential U, the mean first passage time r to transition from a stable equilibrium x to a 


saddle z is (Bovier et al. 2004): 


r = 


2 7 r /1det V 2 17(z)| /2 (U(z) — U(x)) 


| Ai (z)| V det V 2 U(x) 


exp 


cr^ 


(1 + O (cr |log (cr)|)) 


(50) 


33 















V 2 I/(z) is the Hessian of the potential at z, and V 2 H(x) is the Hessian of the potential at 


Ai(z) is the negative eigenvalue of the Hessian at the saddle. Bouchet and Reygner 


(2015) obtained a similar expression has for a non-gradient system with quasi-potential V. 
Their estimate for r is: 


t = 


|detV 2 R(z)| 


2 7T 


|Ai(z)|V detV 2 H(x 


- exp 


2 (V(z)-V(x)) 


O' 


exp 


F(p(t)) dt 


(51) 


Ai(z) is the unstable eigenvalue of the full deterministic skeleton (not just the 
quasi-potential) at the saddle. p{t) is the least action path from x to z: 


p l (t) = VV(p(t)) + Q(p(t)) 


(52) 


F is the divergence of the circulatory component, F(x) = V • Q(x). 

For example 2, we can examine how the mean first passage time estimates from the 


quasi-potential correspond to simulation results (Figure S2). For this example, numerical 
integration along the least action path suggests that J^° F(p(t)) dt ~ 0, so we drop this 
term in the approximation. Note that the quasi-potential approximations closely match the 
means of the simulated first passage times. Of course, there will always be outliers, and 
these can be seen in the tails of the distributions of simulated first passage times in 
figure [S2j 


I Further example 

In the following, we examine a bistable model that illustrates the different ways that 
the stability metrics described in this paper can classify stable states. The model is: 


dX = /, (X, Y) dt + a d\\\ 
dY = f 2 (X, Y) dt + a dW\. 


(53) 


The deterministic skeleton is given by: 
fi(x, y) = -2 ah x exp (- (b\ x 2 + b 2 y 2 )) 
f 2 (x,y) = -2 a b 2 x exp (- (bi x 2 + b 2 y 2 )) 


2d\ (x — c) exp ((di (x — c) 2 + d 2 (y — c) 2 )) 
2 d 2 (y — c) exp ((di (x - c) 2 + d 2 (y - c) 2 )) 


(54) 

This model does not represent any particular ecological process and was instead chosen for 
its ability to illustrate the range of relationships that are possible between the stability 
metrics we discuss. This is a gradient system, with potential function 

U(x,y) = 1 - a exp (- (bix 2 + b 2 y 2 )) - exp (- (d\ (x - c) 2 + d 2 (y — c) 2 )) . (55) 

For all of the parameter values we consider, the system will have two stable states, ei and 
e 2 , separated by a saddle e s . This example will show that each equilibria can be classified 
as more stable by any combination of the stability metrics. Without loss of generality, the 
stable-state with larger x -value, e 2 , will be more stable according to metric 3 (the basin 
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depth metric). 

In case 1, the parameter values are a = 0.9, b\ = 1, b 2 = 1, c = 1.2, d\ = 1.2, d 2 = 1-2. 
e 2 is more stable by all three metrics (Figure S3). 

Metric 1: Re (A (e x )) = -1.16769, Re (A (e 2 )) = -1.75728 
Metric 2: ||ei - ej = 0.643635, ||e 2 - ej = 0.859567 (56) 

Metric 3: U (e s ) - U (ei) = 0.0842552, U (e. s ) - U (e 2 ) = 0.204706 

In case 2, the parameter values are a = 0.9, bi = 2, b 2 = 2, c = 1.8, d\ = 0.8, d 2 = 0.8. 
e 2 is more stable by metric 3 (depth) and metric 2 (basin width) but not metric 1 (linear 
stability) (Figure S4). 

Metric 1: Re (A (e x )) = -3.51331, Re (A (e 2 )) = -1.59979 

Metric 2: ||ei - e s || = 1.05186, ||e 2 - ej = 1.48722 (57) 

Metric 3: U (e s ) - U (e x ) = 0.639468, U (e. s ) - U (e 2 ) = 0.73379 

In case 3, parameter values are a = 0.9, 61 = 1, b 2 = 1, c = 2.5, d\ = 1.2, d 2 = 1.2. e 2 is 
more stable by metric 3 (depth) and metric 1 (linear stability), but not metric 2 (basin 
width) (Figure S5). 

Metric 1: Re (A (ei)) = -1.79998, Re (A (e 2 )) = -2.39984 

Metric 2: ||e x - ej = 1.81859, ||e 2 - ej = 1.71693 (58) 

Metric 3: U (e s ) - U (e x ) = 0.837959, U (e s ) - U (e 2 ) = 0.937962 

In case 4, the parameter values are a = 0.9, b\ = 0.9, b 2 = 0.9, c = 2.37, d\ = 0.78, 
d 2 = 1.46. e 2 is more stable by metric 3 (depth), but not metric 1 (linear stability) or 
metric 2 (basin width) (Figure S 6 ). 

Metric 1 : Re(A(ei)) = —1.61995, Re (A (e 2 )) = —1.55978 

Metric 2: ||ei — e s || = 1.80573, ||e 2 — e s || = 1.77222 (59) 

Metric 3: U (e s ) — U (e^ = 0.81102, U (e s ) — U (e 2 ) = 0.91103 

These four cases show that an equilibrium in a bistable system can be classified as 
“more stable” by any combination of the three metrics. Hence it is important to recognize 
that each metric conveys a different piece of information about stability. 
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(a) CT = 0.02 


(b) a = 0.1 



Figure SI. Approximations to the steady-state probability density of equations Q, obtained 
from a long-time (f = 25000) simulation, with four different noise intensities. Integration was 
performed with the Euler-Maruyama method and At = 0.025. Variables are scaled, so the units are 
dimensionless. The horizontal axis is resource population density and the vertical axis is consumer 
population density. White corresponds to high probability density. The information conveyed in 
each plot depends on the noise intensity (see appendix [Pj). 
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Figure S2. Simulation results for first passage times in example 2. The initial point is eA, and 
the time step is 0.05. 500 realizations were generated at each noise level. The width of each gray 
shape corresponds to the frequency with which each first passage time was observed. The black line 
is the small-noise approximation of the mean first passage time from the formula in appendix [H) 
Note that the small-noise approximation matches the means of the distributions well. At all noise 
levels, the simulations included outliers that escaped from the basin of attraction much faster than 
the small-noise prediction. 
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Figure S3, (a) The potential function for case 1 of the system in appendix [IJ (b) A comparison 

of three different metrics of stability for the system in case 1. 
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Figure S4. (a) The potential function for case 2 of the system in appendix |TJ (b) A comparison 

of three different metrics of stability for the system in case 2. 
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Figure S5. (a) The potential function for case 3 of the system in appendix |TJ (b) A comparison 

of three different metrics of stability for the system in case 3. 
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Figure S6. (a) The potential function for case 4 of the system in appendix [IJ (b) A comparison 

of three different metrics of stability for the system in case 4. 
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